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ABSTRACT 


Large-eddy  simulation  (LES)  output  was  used  to  study  the  terms  of  the  resolved  scale 
turbulent  kinetic  energy  budget  (TKE).  The  computation  produces  a  three-dimensional, 
time  dependent  simulation  of  the  airflow  above  and  within  a  plant  canopy  where  the 
lowest  third  of  the  domain  was  occupied  by  a  drag  layer  and  heat  sources  that  represented 
a  forest.  Shear  and  buoyant  production  of  TKE  computed  from  the  simulated  resolved 
scale  flow  fields  were  principal  sources  in  the  upper  canopy;  the  destruction  of  TKE  due 
to  canopy  drag  effects  and  the  transfer  to  subgrid  scales  (dissipation)  occurred  primarily  in 
the  upper  half  of  the  forest  where  the  foliage  density  was  large;  and  turbulent  transport 
showed  a  loss  at  the  canopy  top  and  a  gain  within  the  canopy.  These  general  features 
have  been  found  in  various  plant  canopies  in  field  experiments,  higher-order  closure 
models,  and  wind  tunnel  studies,  but  in  all  such  studies,  there  is  a  lack  of  information 
concerning  the  pressure  transport  term. 

Previous  studies,  both  experimental  and  numerical,  typically  incorporated  pressure 
effects  into  a  residual  term  due  to  the  difficulty  in  directly  measuring  turbulent  pressure 
fluctuations  and/or  the  lack  of  understanding  on  how  to  parameterize  such  terms.  In  the 
present  LES  study,  the  pressure  was  calculated  directly;  thus,  the  pressure  transport  term 
could  also  be  calculated.  Above  the  canopy,  pressure  transport  appeared  to  balance 
approximately  one-third  of  the  turbulent  transport,  while  near  the  canopy  top  and  below 
pressure  transport  was  the  same  sign  as  turbulent  transport  showing  a  sink  near  the 
canopy  top  and  a  source  below.  The  transport  terms  accounted  for  over  half  of  the  TKE 
sink  at  the  canopy  top,  and  in  the  lowest  two-thirds  of  the  canopy  the  transport  terms 


iii 


were  the  dominant  source  terms  in  the  budget.  Moreover,  the  pressure  transport  was  the 
largest  source  of  turbulent  kinetic  energy  in  the  lowest  levels  of  the  canopy.  These  results 
indicate  that  pressure  transport  is  important  in  the  plant  canopy  turbulent  kinetic  energy 
budget,  especially  in  the  lowest  portion  of  the  stand. 
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1.  Introduction 


Turbulent  kinetic  energy  budgets  are  used  to  study  the  relative  importance  of  different 
physical  processes  that  govern  turbulent  motions.  The  presence  of  a  plant  canopy  in  a 
boundary  layer  provides  unique  challenges  for  measurements,  modeling,  and 
understanding  of  the  physical  processes.  Most  notably,  the  plant  canopy  imposes  an 
aerodynamic  drag  on  the  flow  and  creates  turbulent  motions  in  the  wakes  of  the  plant 
elements.  This  latter  process  is  an  additional  source  of  turbulent  kinetic  energy,  called 
wake  production  (Wilson  and  Shaw,  1977;  Raupach  and  Shaw,  1982;  Raupach  et  al., 
1986),  but  consideration  must  be  given  to  the  scale  of  motion  generated  in  such  wakes.  In 
addition,  the  canopy  may  act  as  a  heat  source  or  sink,  or  both  depending  on  the 
temperature  difference  between  the  layers  of  the  canopy  and  the  surrounding  air  at  each 
level. 

Lesnik  (1974)  provides  an  early  picture  of  the  plant  canopy  turbulent  kinetic  energy 
budget.  Data  from  a  twenty  year  old  pine  forest  were  analyzed  under  different  stability 
conditions.  In  all  three  stability  cases  presented,  the  shear  production  was  a  large  source 
at  the  canopy  top  where  the  wind  shear  gradient  was  the  largest;  the  dissipation  was  a 
large  loss  in  the  upper  canopy,  while  the  turbulent  transport  showed  a  significant  loss  at 
the  top  of  the  canopy  and  a  gain  within  the  canopy.  The  buoyancy  was  small  in 
comparison  with  the  other  terms  and  was  not  shown.  The  unique  feature  of  the  budget 
presented  by  Lesnik  was  the  role  of  the  turbulent  transport  term,  where  the  large  gain  in 
the  lower  canopy  levels  represents  an  import  of  turbulent  energy  fi-om  the  primary 
production  region  in  the  upper  canopy  levels. 
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Analyses  of  recent  field  experiments  (Shi  etal,  1987;  Leclerc  et  al,  1990;  Meyers  and 
Baldocchi,  1991)  support  these  results.  Meyers  and  Baldocchi  (1991)  presented  TKE 
budgets  for  a  deciduous  forest  under  neutral  stability  conditions.  In  the  TKE  budget  the 
measured  profiles  of  shear  production  and  turbulent  transport,  and  the  dissipation,  found 
as  a  residual,  were  in  qualitative  agreement  with  Lesnik  (1974).  Meyers  and  Baldocchi 
(1991)  found  that  wake  production  slightly  exceeds  shear  production  at  all  levels  within 
the  canopy  except  near  the  canopy  top. 

The  effects  of  atmospheric  stability  on  the  TKE  budget  were  investigated  in  Shi  et  al. 
(1987)  and  Leclerc  et  al.  (1990)  using  data  from  a  deciduous  forest  collected  during  the 
Camp  Borden  experiment  (Shaw  et  al.,  1988).  The  most  dramatic  changes  were  observed 
for  stable  conditions  (/j/L>0.25  )  which  were  not  simulated  in  the  present  LES  study.  In 
general,  their  measurements  showed  that  normalized  shear  production  and  dissipation, 
found  as  a  residual,  increase  with  the  onset  of  stability,  and  that  buoyancy  was  a  gain 
under  unstable  stratification  and  a  loss  in  stable  conditions.  However,  their  measurements 
showed  relatively  little  variation  in  the  budget  terms  over  the  range  of  stability  represented 
in  the  present  study. 

Leclerc  et  al.  (1990)  also  presented  the  TKE  budget  in  an  alternative  manner  where 
the  turbulent  kinetic  energy  was  split  into  two  different  scales,  turbulent  shear  kinetic 
energy  (SKE)  and  turbulent  wake  kinetic  energy  (WKE),  following  Shaw  and  Seginer 
(1985).  This  mathematical  split  was  based  on  the  constraint  that  total  TKE  was 
conserved.  In  presenting  the  budgets  in  this  manner,  the  canopy  now  acts  as  a  loss  of 
kinetic  energy  on  the  “large”  shear  producing  scales,  and  a  gain  in  the  “small”  wake  scales. 
This  process  represents  the  action  of  the  drag  elements  to  suppress  turbulent  motions  as 


well  as  the  mean  flow  and,  at  the  same  time,  to  create  small  scale  motions  in  the  wake  of 
plant  elements.  Profiles  showed  that  the  canopy  drag  was  the  primary  destruction  term  in 
the  SKE  budget  while  shear  production  was  the  primary  source  in  the  upper  canopy.  The 
turbulent  transport  was  a  loss  near  the  canopy  top  and  a  gain  within  the  canopy.  Most 
importantly,  the  residual  term  represented  viscous  dissipation,  pressure  transport  and 
accumulated  errors.  Under  unstable  conditions,  the  residual  was  a  gain  within  the  canopy 
and  a  loss  at  the  canopy  top  and  above.  Thus,  pressure  transport  may  be  a  source  of  SKE 
in  the  canopy  and  a  loss  above. 

These  field  experiments  contribute  to  our  understanding  of  turbulent  flow  within  a 
canopy;  however,  the  effects  of  turbulent  pressure  fluctuations  are  still  largely  uncertain. 
Zhuang  and  Amiro  (1994)  provided  stronger  evidence  as  to  the  roles  of  the  pressure 
perturbations  in  plant  canopy  flow.  They  evaluated  the  TKE  budget  during  coherent 
motions  in  a  deciduous  forest  using  data  collected  during  the  Camp  Borden  experiment. 
The  pressure  perturbations  were  not  measured  due  to  the  inherent  difficulty  in  measuring 
static  pressure  perturbations;  rather,  they  were  calculated  by  taking  the  divergence  of  the 
momentum  equation  to  form  a  Poisson  equation  for  the  perturbation  pressure.  The 
pressure  derived  in  this  manner  neglects  the  contributions  from  the  lateral  derivatives  due 
to  the  lack  of  measurements.  Zhuang  and  Amiro  cite  several  studies  that  support  the 
contention  that  for  coherent  motions  lateral  effects  are  small  in  comparison  to  motions  in 
the  x-z  plane,  and  that  pressure  perturbations  derived  from  the  two  dimensional  field 
should  account  for  most  of  the  pressure  effects  in  the  budgets.  In  addition,  the  TKE 
budget  presented  was  for  an  ensemble  average  of  two  dimensional  coherent  structures,  as 


a  result  the  budget  contained  both  longitudinal  and  vertical  pressure  transport  terms. 
Budget  profiles  showed  that  the  vertical  pressure  transport  was  a  loss  at  the  canopy  top 
and  a  gain  within  the  canopy,  and  longitudinal  pressure  transport  was  always  a  loss. 
Although  the  pressure  perturbations  have  been  difficult  to  ascertain  within  a  plant  canopy, 
several  studies  have  measured  surface  pressure  fluctuations. 

Maitani  and  Seo  (1985)  measured  surface  pressure  fluctuations  in  a  wheat  canopy  and 
vertical  velocity  fluctuations  within  and  above  the  canopy.  They  assumed  that  the 
pressure  fluctuations  in  the  canopy  could  be  approximated  by  the  surface  pressure. 
Vertical  plots  of  the  measured  pressure-velocity  covariance  and  an  integration  of  the 
pressure-velocity  cospectrum  over  the  analyzed  frequencies  were,  in  general,  negative 
indicating  that  the  pressure  driven  flux  of  TKE  was  downward.  Their  measurements 
showed  that  the  pressure  transport  term  in  the  TKE  budget  was  not  negligibly  small. 
Shaw  et  al.  (1990)  and  Shaw  and  Zhang  (1992)  used  measured  surface  pressure 
fluctuations  and  velocity  fluctuations  measured  during  coherent  motions  in  a  forest  to 
determine  the  role  of  pressure  fluctuations  in  canopy  turbulence.  Shaw  et  al  (1990) 
compared  surface  pressure  measurements  to  pressure  fluctuations  calculated  through  a 
Poisson  equation  to  demonstrate  that  pressure  fluctuations  at  the  surface  were  primarily 
created  by  velocity  perturbations  in  the  high  shear  region  near  the  top  of  the  forest.  Shaw 
and  Zhang  (1992)  found  that  longitudinal  velocities  measured  in  the  trunk  space  are 
strongly  correlated  with  surface  pressure,  and  that  peak  correlations  occurred  with  near¬ 
zero  time  lag.  Shaw  and  Zhang  stated  that  this  was  strong  evidence,  but  not  proof,  that 
turbulence  at  low  levels  in  the  forest  is  driven  by  pressure  fluctuations.  Pressure  induced 


motions  in  the  lower  canopy  would  account  for  a  significant  pressure  transport  term  in  the 
TKE  budget. 

An  alternative  to  field  experiments  has  been  the  use  of  closure  models  to  investigate 
canopy  budgets  (Wilson  and  Shaw,  1977;  Meyers  and  Paw  U,  1986;  Meyers  and  Paw  U, 
1987;  Meyers  and  Baldocchi,  1991).  Although  these  models  give  qualitative  agreement 
with  the  profiles  of  shear  production,  wake  production,  destruction  effects,  and  turbulent 
transport,  pressure  effects  were  either  assumed  negligible  or  parameterized.  As  a  result, 
no  conclusions  can  be  drawn  regarding  the  role  of  pressure  perturbations  in  the  budgets. 

Wilson  (1988)  presented  a  higher  order  closure  model  using  a  different  formulation. 
The  turbulent  kinetic  energy  budget  was  split  into  two  different  wavebands  in  a  manner 
that  was  similar  to  the  discussion  in  Shaw  and  Seginer  (1985).  The  computed  profiles 
where  in  qualitative  agreement  with  Shaw  and  Seginer.  As  with  other  studies,  the  pressure 
transport  term  was  assumed  negligible,  and  no  conclusions  were  made  on  the  role  of 
pressure  fluctuations  in  the  SKE  budget.  Wilson  did  find  that  the  calculated  level  of  wake 
kinetic  energy  had  no  feedback  on  the  levels  of  shear  kinetic  energy,  mean  velocity,  or 
Reynolds  stress;  therefore,  he  proposed  that  it  was  not  necessary  to  include  wake  kinetic 
energy  in  Ws  closure  model. 

Others  have  turned  to  wind  tunnel  studies  to  investigate  the  airflow  within  and  above  a 
plant  canopy.  Most  notable  are  the  studies  of  Raupach  et  al.  (1986)  and  Brunet  et  al. 
(1994)  for  airflow  in  an  artificial  wheat  canopy.  Budgets  were  presented  that  are  in 
qualitative  agreement  with  field  experiments  and  mathematical  models.  Moreover,  both 
studies  offer  insight  to  the  possible  role  of  the  turbulent  pressure  perturbations  within  and 
above  a  plant  canopy.  In  the  TKE  budget,  Raupach  et  al.  (1986)  found  significant 
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differences  between  two  calculations  of  the  dissipation.  They  suggested  that  the 
difference  may  be  due  to  the  vertical  pressure  transport  indicating  that  pressure  transport 
was  a  source  above  the  canopy  in  approximate  balance  with  turbulent  transport,  and  that 
pressure  and  turbulent  transport  were  both  sources  of  TKE  within  the  canopy.  The  study 
of  Brunet  et  al.  (1994)  also  suggested  that  pressure  transport  was  a  source  above  the 
canopy.  However,  their  residual  analysis  indicated  that  pressure  transport  was  a  loss 
within  the  canopy  which  contradicts  the  earlier  results  of  Raupach  et  al.  (1986).  Brunet  et 
al.  (1994)  cautioned  that  their  analysis  cannot  be  taken  as  proof,  and  due  to  the 
contradictions  with  Raupach  et  al.  (1986)  stated  that  “there  is  clearly  an  urgent  need  for 
high-quality  measurements  of  fluctuating  pressure  in  the  context  of  plant  canopy  flows.” 

LES  may  have  a  partial  answer  to  this  call  for  measurements.  For  example,  Moeng 
(1984)  suggested  that  large-eddy  simulations  (LES)  may  be  a  way  of  obtaining  “data”  to 
study  and  understand  the  physical  processes  within  the  turbulent  planetary  boundary  layer. 
The  LES  described  in  Moeng  (1984)  and  later  modified  in  Moeng  and  Wyngaard  (1988) 
has  been  used  for  several  planetary  boundary  layer  flows  (Moeng  and  Sullivan,  1994). 
The  model  has  been  tested  for  consistency  against  three  LES  codes  written  for  the 
convective  atmospheric  boundary  layer  (Nieuwstadt  et  al.,  1993),  in  which  the  models 
were  run  for  the  same  case  and  compared  favorably.  However,  in  the  TKE  budget  for  the 
simulated  case,  the  pressure  transport  term  was  small  and  the  scatter  between  models 
prevented  any  firm  conclusions,  but  it  was  noted  that  some  of  the  scatter  was  probably 
caused  by  an  insufficient  averaging  period.  Favorable  results  were  obtained  in  Moeng  and 
Wyngaard  (1989)  where  they  compared  their  simulated  TKE  budget  to  budget  terms  fi'om 
a  tank  experiment  (Deardorff  and  Willis,  1985)  and  aircraft  observations  (Lenschow  et  al.. 


1980).  Thus,  LES  proves  to  be  a  valid  technique  to  study  turbulent  flows  within  the 
atmospheric  boundary  layer,  and  it  is  possible  that  LES  may  be  adapted  to  plant  canopy 
airflows. 

Shaw  and  Schumann  (1992)  were  the  first  to  employ  the  LES  technique  to  the  plant 
canopy  environment.  Simulated  mean  velocity  profiles,  vertical  profiles  of  Reynolds 
stress,  turbulent  kinetic  energy,  and  velocity  skewness  qualitatively  matched  observations. 
Recently,  Kanda  and  Hino  (1994)  used  LES  to  explore  coherent  motions  within  and  above 
a  plant  canopy;  however,  the  simulation  was  for  the  developing  stage  of  turbulence  and 
not  fully  developed  turbulent  flow. 

In  this  present  study,  the  LES  of  Moeng  (1984)  and  Moeng  and  Wyngaard  (1988),  as 
modified  by  Patton  et  al.  (1994)  to  include  a  plant  canopy,  was  used  to  generate  model 
data.  The  motivation  was  to  use  the  three  dimensional  output  of  velocity,  temperature, 
subgrid  scale  kinetic  energy,  and  pressure  calculated  through  a  Poisson  equation  to 
calculate  all  terms  in  the  resolved  scale  turbulent  kinetic  energy  budget  with  particular 
emphasis  on  determining  the  role  of  the  pressure  perturbations  in  the  budget.  Several 
simulations  were  run  with  different  canopy  specifications.  Comparisons  were  made 
between  simulations  and  with  field  experiments,  closure  models,  and  wind  tunnel  results. 


2.  Methods 


2.1  Equations 

LES  explicitly  calculates  the  large,  or  resolved,  scales  of  a  turbulent  flow  while  the 
small  scales,  the  subgrid  scales,  are  parameterized.  The  assumption  is  made  that  the 
resolved  scales  contain  most  of  the  energy  and  are  fairly  insensitive  to  the  effects  of  the 
subgrid  scale  parameterizations.  The  resolved  scale  is  mathematically  defined  by  applying 
a  grid  volume  average  or  filter  process  to  the  variables  of  the  governing  equations 
(Leonard,  1974).  In  the  current  model  the  resolved  scale  was  defined  using  a  wave-cut  off 
filter  in  the  horizontal  and  a  grid  volume  average  in  the  vertical  (Moeng  and  Wyngaard, 
1988).  The  derivation  of  the  resolved  scale  equations  was  presented  in  Moeng  (1984); 
however,  some  details  will  be  presented  here  for  completeness,  and  new  terms  will  be 
introduced  to  represent  the  presence  of  the  plant  canopy.  Specifically,  a  form  drag  term 
and  a  canopy  heat  source  are  added  to  the  resolved  scale  momentum  and  energy 
equations,  respectively.  With  the  Boussinesq  approximations,  the  resolved  scale 
conservation  equations  of  mass,  momentum,  and  energy  for  plant  canopy  flow  are 
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where  the  overbar  represents  the  filtering  process.  When  filtering  the  advection  terms,  it  is 
necessary  to  partition  the  flow  variables  into  a  resolved  scale  part  and  a  subgrid  scale  part; 


such  as,  M,.  =Ui+  u'i ,  and,  as  a  result,  additional  terms  arise  in  the  resolved  scale 
equations  (Leonard,  1974;  Moeng,  1984).  This  is  similar  to  additional  terms  arising  when 
forming  the  traditional  Re3molds  decomposition. 
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In  the  momentum  equation  the  additional  terms  appear  as  the  subgrid  scale  stress 
tensor  that  is  written  as 

(4) 

where 

(5) 

In  equation  (4)  the  normal  components  have  been  subtracted  from  the  stress  tensor,  and 
an  equal  and  opposite  term  is  added  to  the  pressure  term  which  has  the  form  of 


(6) 


where  p"  is  the  pressure  deviation  from  the  horizontal  mean  (represented  by  the  double 
prime)  since  the  mean  pressure  field  has  been  separated  from  the  resolved  scale  pressure  in 
equation  (2).  The  third  term  on  the  right  hand  side  of  equation  (6)  is  the  kinetic  energy 
component  of  the  advection  term.  The  substitution  of  (6)  into  the  resolved  scale 
momentum  equation  gives  (see  Appendbc  A  for  more  details  on  the  derivations) 


where  e'^R^jl  is  the  subgrid  scale  (SGS)  kinetic  energy,  and  the  variables  in  the 
conservation  equations  have  the  following  meaning;  is  the  velocity  component  in  the 
Xi  direction,  6  is  the  potential  temperature,  l/^o  is  the  volumetric  expansion  coefficient, 
and  g  is  the  gravitational  acceleration.  The  density  is  given  as  p,  while  (p)  is  the 
horizontal  mean  pressure,  is  the  SGS  shear  stress,  and  is  the  SGS  heat  flux.  The 
canopy  effects  of  form  drag  and  heat  source  are  represented  as  iF)  and  S,  respectively. 


The  resolved  scale  conservation  equations  can  be  solved  with  a  knowledge  of  the  SGS 
kinetic  energy,  the  SGS  fluxes,  and  the  canopy  effects. 

The  subgrid  scale  kinetic  energy  was  solved  using  a  budget  equation  following  the 
formulation  in  Deardorff  (1980).  The  equation  has  the  form 


ae'  -de'  -TT^u.  g-^  du'^e'  +  fjp) 
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where  the  terms  on  the  right  hand  side  represent:  advection  of  SGS  kinetic  energy  by  the 
resolved  scale  velocity,  shear  production  due  to  the  action  of  SGS  Reynolds  stress  on  the 
resolved  scale  velocity  gradient,  SGS  buoyancy  production,  SGS  transport  effects,  and 
dissipation.  Closure  assumptions  are  necessary  to  solve  equation  (8).  Specifically, 


and  heat 


downgradient  diffusion  was  assumed  for  the  SGS  momentum  and  heat 

=  w'O’^  fluxes,  and  for  the  transport  terms.  Kolmogorov’s  hypothesis  was  used  for 


the  dissipation  rate.  It  is  believed  that  the  resolved  scale  motions  are  insensitive  to  these 
parameterizations.  The  last  term  in  (8)  represents  an  enhanced  dissipation  rate  due  to  the 
presence  of  the  plant  canopy  modeled  after  Shaw  and  Schumann  (1992).  The  term  was 
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chosen  to  represent  the  removal  of  SGS  kinetic  energy  by  the  action  of  canopy  drag,  and 
was  based  on  the  assumption  that  eddies  created  in  the  wakes  of  the  plant  elements  are 
much  smaller  than  the  scale  of  the  SGS  eddies;  therefore,  the  wake  eddies  dissipate  rapidly 
making  no  contribution  to  turbulent  kinetic  energy. 

Next,  the  SGS  fluxes  in  the  resolved  scale  equations  were  parameterized  assuming  that 
the  fluxes  can  be  represented  by  the  resolved  scale  strain  rates  in  the  following  manner 
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(10) 


where  the  SGS  eddy  diffusivities,  Kj^  and  ,  were  assumed  functions  of  the  subgrid 

scale  kinetic  energy  and  a  dissipation  length.  Further  details  of  the  subgrid  scale  energy 
budget  and  the  subgrid  scale  fluxes  are  found  in  Moeng  (1984). 

Finally,  canopy  effects,  form  drag  and  heating  source,  must  be  parameterized.  The 
form  drag  of  the  momentum  equation  is  not  simply  a  body  force  that  is  added  to  the 
equation,  rather  the  term  represents  the  combination  of  viscous  drag  and  a  discontinuity  in 
pressure  across  drag  elements  within  the  grid  volume,  and  it  is  analogous  to  the  form  drag 
presented  in  Wilson  and  Shaw  (1977)  and  Raupach  and  Shaw  (1982).  The  form  drag  term 
was  parameterized  in  a  conventional  manner  as  a  product  of  a  drag  coefficient,  cross- 
sectional  body  area,  and  the  square  of  the  time  dependent  velocity.  The  drag  force  in  the 
Xt  direction  is 


F  =  -C,a{z)VR 


(11) 


where  V  is  the  scalar  wind  speed;  Q  is  a  drag  coefficient  assumed  constant  and  equal  to 
0.15  according  to  measurements  in  a  deciduous  forest  (Shaw  et  al,  1988);  and  a{z) 
represents  the  plant  area  density  that  varies  in  the  vertical. 

Since  no  attempt  was  made  to  simulate  the  thermal  energy  budget  of  the  canopy,  the 
canopy  heat  source  term,  represented  by  S,  was  constant  in  time  and  modeled  after  Shaw 
and  Schumann  (1992).  The  magnitude  of  the  canopy  heat  source  was  fixed  at  the  canopy 
top  using  a  specified  canopy  top  heat  flux,  and  the  source  decays  exponentially  as  a 
function  of  an  extinction  coefficient  and  the  downward  cumulative  leaf  area  index.  This 
describes  the  canopy  as  being  the  greatest  source  of  heat  in  the  upper  levels  where  the 
solar  radiation  load  is  the  greatest,  and  as  a  weaker  source  of  heat  in  the  lower  canopy 
where  solar  radiation  is  attenuated. 

2.2  Numerical  model 

The  governing  equations  were  integrated  in  time  using  the  second-order  Adams- 
Bashforth  scheme  with  the  horizontal  derivatives  evaluated  using  a  pseudospectral  method 
(Fox  and  Orszag,  1973),  and  the  vertical  derivatives  evaluated  using  second-order 
centered  finite  differences.  The  boundary  conditions  for  the  conservation  equations  were 
periodic  in  the  horizontal;  the  upper  boundary  was  specified  as  a  fnctionless  rigid  lid  with 
zero  mass,  momentum,  heat,  and  SGS  kinetic  energy  flux;  and  the  bottom  boundary 
employed  the  no-slip  condition.  The  time  step  for  the  simulations  was  equal  to  0.1 
seconds  and  the  grid  spacing  of  2  meters  was  equidistant  in  all  directions.  There  were  96 
X  96  horizontal  grids  and  31  vertical  grids  (representing  a  domain  size  of  192m  x  192m  x 
60m)  with  the  lowest  one-third  of  the  domain  (10  grid  points)  occupied  by  a  20m  tall 


forest.  The  forest  was  defined  as  a  vertical  distribution  of  leaf  area  density.  Figure  1  gives 
the  two  profiles  used  in  this  study.  The  profiles  are  defined  as  a  sparse  canopy,  leaf  area 
index  (LAI)  equal  to  2,  and  a  dense  canopy  with  LAI  =  5.  The  leaf  area  index  is  the 
vertical  integration  of  the  leaf  area  density  from  the  surface  to  the  canopy  top. 

The  conservation  equations  were  solved  on  a  non-staggered  horizontal  grid  and  a 
staggered  vertical  grid.  The  vertical  velocity,  w ,  the  subgrid  scale  kinetic  energy,  e' ,  and 
the  leaf  area  density  are  defined  at  the  surface  and  on  equally  spaced  grids  above  (vv 
level).  The  streamwise  velocity,  u  ,  the  lateral  velocity,  v,  potential  temperature,  0 ,  and 
the  dynamic  pressure,  P* ,  are  defined  at  a  half  grid  point  above  the  surface  and  the 
intermediate  levels  above  (m  level).  The  budget  terms  were  calculated  at  the  u  level 
resulting  in  the  budget  terms  being  offset  from  the  location  of  the  leaf  area  density. 

The  LES  was  run  for  several  different  cases  where  the  simulated  flow  was  forced  by 
treating  the  horizontal  mean  pressure  gradient  as  an  external  forcing.  The  pressure 
gradient  was  adjusted  at  each  time  step  to  ensure  a  constant  mass  flux  across  the  upwind 
y-z  plane  where  the  mass  flux  was  determined  by  defining  a  constant  mean  wind  speed, 
U ,  at  the  upwind  boundary.  The  mean  wind  speed  was  varied  between  simulations  in 
order  to  study  the  environmental  influences  on  the  budget.  The  canopy  density  and  the 
strength  of  the  heat  sources  were  also  varied  between  simulations  to  determine  canopy 
influences.  As  mentioned  previously,  two  canopy  densities  were  used  in  this  study.  In 
addition,  the  canopy  heating  source  was  represented  as  either  a  high  sensible  heat 
exchange  or  a  low  sensible  heat  exchange  by  specifying  two  different  canopy  top  heat 
fluxes,  0,  =  0.125  and  0.005  mKs~\  Table  1  summarizes  environmental  and  canopy 


forcing  used  in  the  simulations.  The  simulations  are  ordered  in  terms  of  a  stability 
parameter  hJL  (Shaw  et  al,  1988;  Leclerc  et  al,  1990)  where  h  is  the  height  of  the 
canopy  and  L  is  the  Monin-Obukov  length  evaluated  at  the  canopy  top. 

After  the  flow  reached  equilibrium  for  a  simulation,  data  sets  were  saved  at  a  specified 
time  interval  (250s  or  500^).  A  data  set  contains  the  three  dimensional  output  of  the 
resolved  scale  velocity  (m,v,vv),  the  resolved  scale  potential  temperature  (0),  the  SGS 

kinetic  energy  (e'  ),  and  the  dynamic  pressure  (P*).  However,  to  form  the  budgets  it 
was  convenient  to  recalculate  the  pressure  perturbation  field  using  a  Poisson  equation  of 
the  form 


fp")  d  'Ide' 


dx  “  e  "  dx.  ' 


(12) 


where  the  Poisson  equation  was  obtained  by  taking  the  divergence  of  the  resolved  scale 
momentum  equation.  At  this  point  all  terms  of  the  resolved  scale  TKE  budget  can  be 
calculated. 


2.3  The  TKE  budget  equation 

The  derivation  of  the  resolved  scale  TKE  budget  was  similar  to  the  traditional  free 
atmosphere  averaging  procedures  (Stull,  1988).  Briefly,  a  horizontal  mean  momentum 
equation  was  derived  by  decomposing  the  flow  variables  in  the  resolved  scale  momentum 
equation  (7)  into  a  horizontal  mean  and  a  horizontal  deviation.  The  horizontal  mean 
momentum  equation  was  obtained  by  taking  the  horizontal  average  of  the  decomposed 


equation  (see  Appendix  B).  A  horizontal  perturbation  momentum  equation  was  then 
formed  by  subtracting  the  horizontal  mean  momentum  equation  from  the  resolved  scale 
momentum  equation  (Appendix  C).  The  canopy  drag  term  was  included  in  this  analysis 
resulting  in  drag  also  reducing  perturbation  momentum.  Thus,  a  TKE  budget  derived 
from  the  perturbation  momentum  equation  also  contains  a  term  that  acts  to  remove 
resolved  scale  TKE. 

Wilson  and  Shaw  (1977)  specifically  caution  against  such  a  formulation.  They  stated 
that  adding  a  body  force  to  represent  canopy  drag  would  satisfy  the  momentum  equation, 
but  such  a  term  would  not  adequately  describe  other  canopy  effects  (wake  production) 
since  drag  would  reduce  turbulence  as  well  as  mean  momentum  from  the  flow.  Instead, 
proper  consideration  of  horizontal  averaging  and  spatial  differentiation  could  describe  the 
canopy  effects  in  the  momentum  equation  and  the  TKE  budget.  The  canopy  effects  would 
appear  as  a  form  drag  term  in  the  momentum  equation,  and  the  wake  production  process 
would  naturally  arise  due  to  mathematical  manipulation  of  the  terms  in  the  TKE  budget 
(Wilson  and  Shaw,  1977;  Raupach  and  Shaw,  1982;  Raupach  etal,  1986). 

In  the  present  study,  the  above  is  not  a  concern  since  the  numerical  simulation 
specifically  separates  large  and  small-scale  motions,  and  motions  created  in  wakes  behind 
canopy  elements  are  assumed  to  reside  in  unresolved  subgrid  scales.  Through  the 
momentum  equation,  the  canopy  drag  term  acts  to  suppress  resolved  scale  motions,  and 
this  influence  is  automatically  carried  through  to  resolved  scale  turbulent  kinetic  energy. 
The  filter  process  defines  the  larger  scales  of  the  flow  and  sets  a  lower  limit  on  the  size  of 
the  resolvable  eddies.  Consequently,  it  was  assumed  that  the  wake  production  process 


should  occur  on  the  subgrid  scales,  and  a  “wake  destruction”  process  should  occur  on  the 
resolved  scales. 

This  assumption  is  similar  to  the  ideas  presented  in  Shaw  and  Seginer  (1985).  Shaw 
and  Seginer  stated  that  the  total  TKE  could  be  viewed  to  consist  of  two  components,  the 
kinetic  energy  of  the  shear  generated  eddies  and  the  kinetic  energy  of  the  eddies 
characteristic  of  the  wakes  of  the  canopy  elements.  This  separation  was  based  on  the 
characteristic  length  scales  of  the  eddies.  The  shear  scale  eddies  had  a  length  scale  about 
one  order  of  magnitude  larger  than  the  height  above  ground,  and  the  wake  eddies 
occurred  on  the  scale  of  the  plant  elements.  Shaw  and  Seginer  state  that  there  can  be  two 
orders  of  magnitude  difference  between  the  scales.  This  separation  resulted  in  the  canopy 
having  opposing  effects  on  the  budgets,  where  the  canopy  acted  as  a  sink  for  shear  kinetic 
energy  and  a  source  for  wake  kinetic  energy.  In  the  present  study,  the  resolved  scales  are 
analogous  to  the  shear  scale  and  the  wake  motions  occur  on  the  subgrid  scales.  However, 
as  stated  above,  the  assumption  is  made  that  wake  production  does  not  contribute 
significantly  to  the  subgrid  scales  because  wake  scales  are  small  and  rapidly  dissipate  to 
heat. 

The  resolved  scale  TKE  budget  written  in  mixed  notation  for  steady  state,  horizontally 


homogeneous  conditions  is 
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where  {E"')  =  {ui  ul'l'i)  is  the  horizontal  averaged  resolved  scale  turbulent  kinetic  energy 

(see  Appendix  D  for  more  details  on  the  derivation  of  the  TKE  budget).  The  angled 
brackets  represent  the  horizontal  average  and  the  double  prime  represents  the  deviation 
therefrom.  The  perturbation  momentum  is  defined  as  u"=  m,-  -(«))•  ^  is  the  mean  shear 

production  representing  a  conversion  of  mean  flow  kinetic  energy  to  the  resolved  scale 
turbulent  kinetic  energy.  This  term  is  usually  a  source  and  is  due  to  the  interaction  of  the 
mean  velocity  gradients  with  the  turbulent  velocity  field,  fj,  acts  as  either  a  buoyant 

production  or  destruction  depending  on  the  sign  of  the  vertical  heat  flux.  The  heat  flux  is 
positive  for  unstable  conditions  representing  a  conversion  of  potential  energy  to  TKE;  and 
is  a  sink  for  TKE  under  stable  conditions.  is  the  resolved  scale  turbulent  transport, 

and  Tp  represents  the  vertical  transport  of  resolved  scale  kinetic  energy  by  the  pressure 

fluctuations.  D^d  is  the  parameterized  influence  of  canopy  drag.  This  term  represents  the 

rate  of  work  performed  by  the  velocity  perturbations  against  the  drag  forces,  and  is 
identical  to  the  canopy  drag  term  in  Shaw  and  Seginer  (1985).  Dggg  represents  the 

combined  effects  of  subgrid  scale  diffusion  of  resolved  scale  TKE  (first  two  terms)  and  the 


transfer  of  resolved  scale  TKE  to  the  subgrid  scales  (last  term).  The  subgrid  scale 
diffusion  effects  were  separated  from  Dgg^  and  added  to  the  resolved  scale  turbulent 

transport  (  7} ),  but,  in  general,  their  contribution  was  small.  To  avoid  confusion  with  the 

resolved  scale  turbulent  transport  ( 2} ),  the  sum  of  resolved  scale  turbulent  transport  and 

subgrid  scale  diffusion  is  referred  to  as  turbulent  transport  ( 7T ).  The  transfer  of  resolved 
scale  TKE  to  the  subgrid  represents  a  conversion  of  energy  to  the  dissipation  scales 
(Moeng,  1984;  Mason,  1994).  In  this  study,  this  is  referred  to  as  subgrid  scale  dissipation. 

The  budgets  presented  here  are  different  from  typical  budget  equations  for  the  canopy 
layer  (Wilson  and  Shaw,  1977;  Raupach  and  Shaw,  1982;  Raupach  etal.,  1986).  First,  as 
noted  above,  there  is  no  wake  production  term  based  on  the  belief  that  this  process  cannot 
be  resolved  in  the  current  LES.  This  assumption  is  based  on  the  fact  that  the  grid 
resolution  is  two  meters  in  the  horizontal  and  vertical;  thus,  only  eddies  of  wavelength 
greater  than  four  meters  can  be  resolved.  Eddies  of  this  size  are  significantly  larger  than 
the  size  of  the  eddies  expected  of  the  wake  production  process  that  occur  on  the  scales  of 
the  plant  elements  (1-50  cm).  Secondly,  there  are  no  dispersive  fluxes  (Raupach  and 
Shaw,  1982;  Raupach  et  al.,  1986).  The  dispersive  fluxes  represent  spatial  correlations 
between  time  averaged  quantities;  for  example,  at  the  edge  of  a  tree  stand  a  consistent 
downdraft  may  occur  on  the  time  average.  This  downdraft  may  be  correlated  with 
turbulent  motions  upwind  of  the  edge,  and  the  net  effect  would  be  a  significant  dispersive 
flux  in  the  TKE  budget.  Dispersive  fluxes  are  absent  from  the  equations  because  the 
budgets  are  for  only  a  spatial  average,  as  opposed  to  performing  a  spatial  or  volume 
average  in  conjunction  with  a  time  average;  however,  an  ensemble  average  of  each  term 


was  formed  by  averaging  the  terms  over  the  number  of  available  data  sets  within  each 
simulation. 


3.  Results  and  Discussion 


3.1  TKE  budget  profile 

The  normalized  budget  for  the  near  neutral  simulation  (S5)  is  given  in  Figure  2.  The 
overall  features  of  the  budget  are  in  agreement  with  previous  studies.  Above  the  canopy 
there  is  an  approximate  local  balance  between  the  shear  production  and  the  rate  of  transfer 
of  energy  to  the  subgrid  scales  (dissipation).  However,  at  the  top  of  the  domain  shear 
production  approaches  zero  due  to  the  specification  of  a  rigid  upper  boundary  with  no 
fluxes  across  the  boundary.  The  subgrid  scale  dissipation  remains  an  expected  loss 
throughout  the  domain  resulting  in  a  local  imbalance  between  production  and  dissipation. 
The  imbalance  is  compensated  for  by  turbulent  transport  that  is  unable  to  export  kinetic 
energy  through  the  upper  lid.  In  previous  experimental  studies,  turbulent  transport  was  a 
maximum  loss  at  the  canopy  top  and  remained  a  loss  above  the  canopy  decreasing  in 
magnitude  away  from  the  top  of  the  canopy.  The  positive  values  of  turbulent  transport 
above  the  canopy  in  this  case  appear  to  be  an  artifact  of  the  boundary  conditions. 

Above  z/^»1.3  the  pressure  transport  is  roughly  opposite  the  turbulent  transport, 

and,  in  regions  where  the  terms  are  not  negligibly  small,  the  magnitude  of  pressure 
transport  is  approximately  30%  of  the  turbulent  transport.  The  relative  percentage  of 
pressure  transport  to  turbulent  transport  varies  somewhat  between  simulations.  This 
feature  cannot  be  considered  conclusive  proof  that  pressure  transport  opposes  turbulent 
transport  above  the  canopy  because  of  the  possible  boundary  influences;  however,  the 
equilibrium  between  production  and  dissipation  with  a  balance  between  the  transport 
terms  was  inferred  in  the  canopy  wind  tunnel  study  of  Raupach  et  al.  (1986). 


There  is  considerable  evidence  that  this  balance  occurs  within  the  unstable  atmospheric 
surface  layer;  for  example,  McBean  and  Elliott  (1975)  measured  turbulent  velocity 
components,  temperature,  and  pressure  at  one  height  in  the  surface  layer.  The  vertical 
profiles  of  the  TKE  budget  terms  were  then  inferred  from  Monin-Obukhov  similarity  using 
the  measured  data.  Their  results  showed  that  the  pressure  transport  was  approximately 
equal  to  and  opposite  the  turbulent  transport.  Brost  et  al.  (1982)  summarized  several 
studies  of  the  surface  layer  and  the  convective  boundary  layer  that  showed  the  pressure 
transport  opposing  turbulent  transport.  In  the  reported  studies,  pressure  transport  was 
determined  as  a  residual  and  the  percentage  of  pressure  to  turbulent  transport  varied  fi-om 
50%  in  Lenschow  et  al.  (1980)  to  100%  in  Wyngaard  and  Cote  (1971),  Rayment  and 
Chaughey  (1977),  and  Chaughey  and  Wyngaard  (1979). 

McBean  and  Elliott  (1975)  proposed  a  simple  model  for  the  observed  pressure 
transport  profile  based  on  measurements  that  showed  the  covariance  (vv " p”) ,  using  the 
current  notation,  and  its  vertical  divergence  were  negative.  These  results  indicated  that 
the  pressure  and  vertical  velocity  perturbations  were  of  opposite  sign  and  that  the  pressure 
transport  was  a  source  of  turbulent  energy  in  the  lower  levels.  McBean  and  Elliott 
postulated  that  a  local  increase  in  perturbation  pressure  occurred  below  the  maximum 
downward  motion  for  a  downward  moving  eddy  in  contact  with  the  boundary,  and  the 
reverse  occurred  for  an  upward  moving  eddy  that  was  in  contact  with  the  boundary.  Both 
motions  result  in  a  negative  covariance,  and  in  order  to  have  a  negative  vertical 
divergence,  the  covariance  must  become  increasingly  negative  with  height.  McBean  and 
Elliott  do  not  comment  on  why  this  should  occur;  however,  Wilczak  and  Businger  (1984) 


found  that  {w"  p")  was  increasingly  negative  with  height  using  tower  measured  data  for  a 
series  of  temperature  ramp  structures. 

Wilczak  and  Businger  (1984)  hypothesized  that  the  pressure  field  retrieved  from  the 
momentum  equations  during  the  ramp  structures  should  account  for  the  pressure 
covariance  term  in  the  turbulent  kinetic  energy  budget,  and  that  the  covariance  became 
more  negative  with  height  was  due  to  observed  phase  differences  between  pressure  and 
vertical  velocity  associated  with  the  measured  inclined  ramp  structures.  As  a  result, 
budget  equations  showed  that  pressure  transport  was  a  source  of  turbulent  kinetic  energy 
in  approximate  balance  with  turbulent  transport.  At  this  point  it  cannot  be  stated  that 
ramp  structures  are  responsible  for  the  observed  pressure  patterns  in  the  studies  other  than 
Wilczak  and  Businger  (1984).  Nevertheless,  pressure  transport  appears  to  balance  a 
percentage  of  the  turbulent  transport  in  the  atmospheric  surface  layer  and  in  the  wind 
tunnel  study  of  Raupach  et  al.  (1986).  It  also  appears  that  such  a  feature  may  occur  above 
the  simulated  canopy  presented  here. 

This  result  has  been  found  in  several  numerical  studies  of  the  atmospheric  boundary 
layer  (e.g.  Deardorff,  1980;  Mason,  1989;  Moeng  and  Wyngaard,  1989).  Deardorff 
(1980)  simulated  the  stratus  capped  convective  mixed  layer  and  found  that  pressure 
transport  offset  one-third  to  one-half  turbulent  transport  except  near  the  cloud  base  where 
the  two  terms  were  both  sources  of  turbulent  kinetic  energy.  Similar  profiles  were  found 
in  the  convective  boundary  layer  simulations  of  Moeng  and  Wyngaard  (1989)  and  Mason 
(1989)  where  pressure  transport  roughly  opposed  turbulent  transport  throughout  the 
boundary  layer  except  just  below  the  inversion  base. 


The  approximate  local  equilibrium  observed  above  the  canopy,  does  not  exist  near  the 
canopy  top  and  below.  In  the  upper  canopy  levels,  shear  production  is  the  dominant 
source  of  turbulent  energy  due  to  the  simulated  large  mean  wind  gradients.  At  the  canopy 
top  shear  production  is  balanced  by  the  destruction  effects  (canopy  drag  effects  and  SGS 
dissipation)  and  the  transport  processes.  In  fact,  just  above  the  canopy  (z/h=l.05) 
turbulent  transport  and  pressure  transport  account  for  over  50%  of  the  TKE  loss.  Table  2 
gives  a  summary  of  the  ratios  between  the  transport  terms  at  the  top  of  the  canopy  to  the 
combined  peak  production.  In  all  simulations,  the  transport  terms  account  for  a  significant 
portion  of  the  TKE  loss. 

Balance  within  the  canopy  occurs  primarily  between  the  loss  due  to  the  canopy  drag 
effects  and  the  gain  due  to  shear  production,  turbulent  transport,  and  pressure  transport. 
However,  shear  production  decreases  quite  rapidly  becoming  negligible  below  z//i=0.65; 
thus,  turbulent  transport  and  pressure  transport  serve  as  the  primary  TKE  sources  in 
roughly  the  lowest  two-thirds  of  the  canopy.  Moreover,  pressure  transport  is  the  primary 
source  below  zlh»05  accounting  for  over  70%  of  the  TKE  source  in  the  lowest  third  of 
the  canopy.  The  significance  of  this  effect  is  difficult  to  visualize  when  the  transport 
terms  are  plotted  on  the  scale  of  the  largest  budget  terms;  therefore,  in  order  to  better 
discern  the  effects  of  the  budget  terms,  the  individual  terms  have  been  plotted  separately  in 
the  following  section. 

3.2  TKE  budget  terms 

Shear  production  is  the  primary  TKE  production  process  in  all  simulations.  Figures  3  a 
and  3b  show  normalized  shear  production  profiles  for  the  simulations  with  leaf  area  index 


of  2  and  5,  respectively.  The  plot  labels  in  the  lower  left  of  each  figure  are  ordered  in 
terms  of  increasingly  negative  values  of  the  stability  parameter  hjL ;  thus,  in  Figure  3  a 

case  SI  is  the  least  unstable,  while  case  S4  is  for  the  most  unstable  simulation.  For  ease  of 
comparison  the  simulations  with  /i/L>-02  will  be  referred  to  as  the  near  neutral 
simulations,  and  when  hjL  <-02  the  simulations  will  be  referred  to  as  the  unstable 
simulations. 

The  general  features  of  the  normalized  shear  production  for  the  sparse  canopy  (LAI  = 
2)  show  a  peak  in  production  at  z/A=0.95  and  a  fairly  rapid  decrease  in  the  canopy 
(Figure  3  a).  The  near  neutral  cases  show  some,  albeit  small,  scatter  among  the 
simulations,  and  below  z//j«0.80  the  profiles  collapse  to  a  similar  profile.  Moreover,  the 
profiles  for  simulations  SI  and  S3  are  in  very  good  agreement  (approximately  5%) 
throughout  the  domain.  These  simulations  were  for  the  same  low  canopy  top  heat  flux, 
and  the  mean  wind  speed  was  varied  from  2  ms~^  in  SI  to  1  ms~^  in  S3.  Clearly  these 
two  profiles  scale  well  when  normalized  in  the  conventional  manner  by  h/ul .  A  slightly 

larger  difference  is  noted  when  comparing  simulations  SI  and  S3  with  simulation  S2 
(approximately  15%).  In  the  latter  case,  the  canopy  top  heat  flux  was  large  and  the  mean 
wind  speed  was  increased  to  4  Although  the  simulated  profiles  are  within  10% 

agreement,  the  profiles  suggest  that  the  peak  production  is  smaller  for  the  high  wind  speed 
case  which  is  physically  unrealistic.  Dimensional  plots  (not  shown)  reveal  the  expected 
order  where  the  shear  production  rate  is  largest  for  the  high  wind  speed  case  and 
decreases  as  a  function  of  decreasing  mean  wind  speed.  The  difference  here  may  be  due 
to  the  number  of  data  sets  analyzed.  Time  series  traces  of  the  canopy  top  fnction  velocity 


showed  considerable  variability  (Patton,  personal  communication),  and,  as  a  result,  the 
number  of  data  sets  may  not  be  sufficient  to  account  for  these  variations.  Small  changes  in 
fiiction  velocity  will  have  significant  effect  on  profiles  since  the  fnction  velocity  is  cubed 
when  normalizing. 

The  largest  differences  in  normalized  production  are  observed  between  the  near 
neutral  simulations  and  the  unstable  simulation  (S4).  The  unstable  case  had  a  high  canopy 
top  heat  flux  with  a  specified  low  mean  wind  speed.  Normalized  shear  production  is 
reduced  by  approximately  40%.  Leclerc  et  al.  (1990)  found  similar  results  where 
normalized  shear  production  increased  with  the  onset  of  stability.  Based  on  the  these 
results,  one  may  conclude  that  shear  production  is  smaller  for  the  unstable  case;  however, 
this  may  not  be  true  for  all  physical  situations.  For  example,  consider  simulations  S3  and 
S4  where  both  simulations  had  a  specified  low  mean  wind  speed  ( 1  w  5"* ),  and  the  only 
difference  between  simulations  was  the  canopy  top  heating.  In  dimensional  form  (not 
shown)  the  shear  production  was  larger  for  the  unstable  simulation  (S4).  This  is  a 
consequence  of  the  enhanced  vertical  motions  due  to  the  larger  canopy  heating. 
Although,  the  enhanced  vertical  motions  will  reduce  the  mean  wind  gradient,  the  vertical 
motions  will  also  increase  the  flux  of  momentum  within  the  canopy  layer,  and  it  appears 
that  this  increase  in  momentum  flux  is  responsible  for  the  larger  shear  production 
(dimensional  form)  observed  in  the  unstable  simulation.  Added  proof  that  the  momentum 
flux  has  increased  is  seen  by  comparing  the  canopy  top  fnction  velocities  (Table  1)  where 
the  fnction  velocity  for  S4  is  approximately  40%  larger  than  simulation  S3. 


The  nondimensional  profiles  reveal  the  relative  importance  of  shear  production  for  a 
given  simulation,  and  the  small  scatter  between  near  neutral  simulations  indicates  that 
similar  physical  processes  account  for  turbulent  kinetic  energy  production.  For  example,  a 
comparison  of  peak  shear  production  (for  both  normalized  and  dimensional  form)  to  total 
peak  TKE  production  reveals  that  shear  accounts  for  over  98%  of  the  total  TKE 
production  in  the  near  neutral  simulations,  while  in  the  unstable  case  shear  production  only 
accounts  for  approximately  75%  of  the  total  peak  production. 

The  normalized  shear  production  profiles  for  the  dense  canopy  (LAI  =  5)  cases  show 
features  similar  to  the  sparse  canopy.  Normalized  peak  production  occurs  at  zlh=0.95 
(Figure  3b),  and  the  near  neutral  profiles  show  some  scatter  collapsing  to  a  single  profile 
near  z//i«0,85 .  In  addition,  the  normalized  peak  production  is  reduced  by  approximately 

30%  for  the  unstable  case,  and,  for  both  simulated  canopies,  the  profiles  are  not  ordered 
strictly  by  the  stability  parameter,  hjL .  Although  the  profiles  for  both  leaf  area  densities 
are  similar,  there  are  several  notable  differences  due  to  the  change  in  canopy  density. 

Normalized  shear  production  is  larger  for  the  dense  canopy,  and,  due  to  the 
aforementioned  variations  in  fiiction  velocity,  a  comparison  was  made  with  dimensional 
plots  (not  shown).  The  dimensional  plots  showed  that  the  shear  production  rate  was 
larger  when  only  considering  a  change  in  canopy  density.  For  example,  the  peak 
dimensional  production  for  simulation  S7  (dense  canopy)  was  approximately  50%  larger 
than  simulation  S3  (sparse  canopy),  and  both  production  rates  (in  dimensional  form)  were 
much  smaller  than  the  other  cases  because  the  canopy  heating  and  the  mean  wind  speed 
were  small.  Further  effects  on  shear  production  due  to  the  change  in  canopy  density  can 


be  seen  when  comparing  the  near  neutral  simulations,  where  the  normalized  peak 
production  value  for  the  sparse  canopy  is  approximately  6.5  while  in  the  dense  canopy  the 
peak  is  about  8.5  representing  a  30%  increase.  A  similar  increase  is  noted  for  the  unstable 
simulations.  These  results  indicate  that  the  mean  wind  gradient  is  larger  for  the  dense 
canopy.  This  is  expected  since  the  dense  canopy  will  extract  more  momentum  due  to  a 
greater  canopy  drag  resulting  in  a  larger  difference  between  the  canopy  top  and  ffee- 
stream  wind  speed;  however,  the  changes  here  may  be  influenced  by  the  specification  of  a 
constant  mass  flux  across  the  upwind  edge  of  the  model  domain.  For  a  given  mean  wind 
speed,  a  constant  mass  flux  across  the  domain  requires  a  stronger  flow  above  the  dense 
canopy  since  the  dense  canopy  will  reduce  the  airflow  within  the  canopy  layer,  and  the 
result  being  a  larger  mean  wind  gradient. 

The  normalized  shear  production  profile  in  the  dense  canopy  also  shows  a  much  more 
rapid  decrease  within  the  canopy.  The  dense  canopy  shear  production  becomes  negligible, 
less  than  10%  of  the  canopy  top  peak,  at  a  slightly  higher  level  {zlh  »  0.65)  than  in  the 
sparse  canopy  (zlh»055).  Meyers  and  Baldocchi  (1991)  also  found  that  normalized 
shear  production  decreased  rapidly  inside  a  deciduous  forest  with  a  leaf  area  index  of  5. 
Their  measurements  showed  production  to  be  negligible  below  z//iw0.75 .  The  difference 

is  likely  due  to  the  canopy  architecture  because  the  foliage  was  concentrated  more  in  the 
upper  levels  of  the  forest  reported  in  Meyers  and  Baldocchi  (1991). 

Figure  4  presents  the  profiles  for  normalized  buoyant  production  where  the  profiles 
are,  as  expected,  ordered  in  terms  of  the  stability  parameter.  In  general  buoyancy  is  a 
negligible  source  when  compared  to  shear  production  except  in  the  unstable  simulations 


S4  and  S8.  Specifically,  in  the  simulated  unstable  cases  the  peak  buoyant  production 
accounts  for  approximately  25%  of  the  total  peak  production,  and  less  than  2%  of  the 
total  peak  production  for  the  near  neutral  simulations.  The  budget  profiles  presented  in 
Leclerc  et  al.  (1990)  showed  similar  results  where  the  normalized  buoyant  production  was 
significant  in  their  unstable  case,  and  normalized  buoyant  production  was  negligible  for 
their  neutral  measurements.  The  common  features  for  both  simulated  canopies  show  that 
buoyant  production  is  a  maximum  at  zlh^  « 1.05 ,  and  is  a  source  at  all  levels  except  for 

cases  S8,  S6,  and  S2  where  buoyancy  is  a  sink  below  z//i«0.45 .  In  the  latter  two  cases 
the  magnitude  of  the  normalized  sink  is  on  the  order  of  10'^. 

The  maximum  buoyant  production  just  above  the  canopy  indicates  that  the  vertical 
heat  flux  was  large  and  positive  at  this  level.  This  result  is  a  direct  consequence  of  the 
specified  canopy  heating  profile.  In  the  simulations,  the  canopy  heating  was  largest  at  the 
top  of  the  canopy  and  the  heating  exponentially  decayed  with  decreasing  height. 
Physically  this  represents  solar  radiation  heating  the  leaves  in  the  upper  canopy  which  in 
turn  heats  the  surrounding  air  given  that  no  transpiration  occurs.  The  decay  into  the 
canopy  represents  the  strong  attenuation  of  solar  radiation  resulting  in  the  canopy  acting 
as  a  weaker  heat  source. 

The  reversal  of  buoyancy  fi'om  a  production  process  to  a  destruction  process  occurs  in 
the  simulations  that  have  a  specified  high  canopy  heat  source  (except  case  S4),  and  is  only 
significant  in  case  S8.  Buoyant  destruction  occurs  under  thermally  stable  conditions  when 
the  vertical  heat  flux  becomes  negative.  A  stable  layer  develops  within  the  simulated 
forest  due  to  the  specification  of  an  elevated  heat  source  resulting  in  the  air  in  the  upper 


canopy  to  be  the  warmest  with  cooler  air  temperature  below.  Such  a  stable  layer  has  been 
found  to  be  typical  of  the  daytime  temperature  profile  within  a  canopy  (Kaimal  and 
Finnigan,  1994);  however,  the  observed  canopy  heat  flux  is  generally  positive  representing 
a  countergradient  flux  of  heat.  The  simulated  downward  heat  flux  may  be  an  inadvertent 
consequence  of  the  boundary  conditions.  The  surface  temperature  was  specified  using  a 
similarity  relationship  and  the  temperature  at  the  lowest  u  level  (a  half  grid  point  from  the 
surface),  while  the  upper  boundary  did  not  allow  a  flux  of  heat  across  the  boundary.  As  a 
result  heat  was  continually  added  to  the  system  and  thermal  equilibrium  was  never 
reached.  The  continual  source  of  heat  in  the  upper  canopy  could  result  in  the  simulated 
downward  heat  flux;  however,  the  simulated  profiles  are  not  unrealistic  since  the 
simulations  are  on  the  order  of  lO’s  of  minutes  and  thermal  equilibrium  may  not  be 
reached  on  this  time  scale  for  certain  atmospheric  conditions. 

The  effect  of  the  mean  wind  shear  on  buoyant  production  is  seen  by  comparing 
simulation  S4  to  simulation  S2  in  the  sparse  canopy  and  by  comparing  simulation  S8  to 
simulation  S6  for  the  dense  canopy.  In  these  simulations  the  canopy  top  heat  flux  was 
large  (0.125  mKs~^),  and  the  mean  wind  speed  was  equal  to  1  in  simulation  S4 
(S8)  with  the  mean  wind  increased  to  4  m  in  simulation  S2  (S6).  The  canopy  heat 
source  is  large,  yet  for  the  high  wind  case  normalized  buoyant  production  is  small. 
Conversely,  buoyant  production  was  large  when  the  canopy  heating  was  large  and  mean 
wind  small.  These  results  may  seem  contradictory  because  for  a  given  leaf  area  index  and 
canopy  heating  rate  one  would  expect  that  the  rate  of  buoyant  production  in  simulation  S4 
(S8)  to  equal  buoyant  production  for  simulation  S2  (S6).  In  dimensional  form  (not 


shown)  the  production  rates  are  approximately  equal,  and,  as  with  normalized  shear 
production,  the  profiles  in  Figure  4  show  the  relative  importance  of  normalized  buoyant 
production  to  total  production  within  each  simulation. 

The  production  processes  are  normally  balanced  by  the  destruction  effects  in  the  TKE 
budget.  In  the  resolved  scale  TKE  budget  there  are  two  destruction  mechanisms,  the  SGS 
dissipation  and  plant  canopy  drag  effects.  The  SGS  dissipation  represents  the  normal 
inertial  cascade  whereby  large  eddies  are  broken  down  into  successively  smaller  eddies 
that  are  ultimately  acted  upon  by  viscosity.  On  the  other  hand,  the  canopy  drag  effects 
represent  the  direct  removal  of  turbulent  energy  due  to  the  interaction  of  resolved  scale 
turbulence  with  the  plant  canopy.  This  process  results  in  a  short  circuit  of  the  normal 
eddy  cascade  (Shaw  and  Seginer,  1985). 

Figure  5a  presents  the  normalized  SGS  dissipation  for  the  sparse  canopy  simulations. 
The  dissipation  is  large  at  the  canopy  top  and  decreases  within  the  canopy.  As  with  shear 
production,  the  dissipation  is  not  strictly  ordered  by  the  stability  parameter.  The  near 
neutral  simulations  show  little  scatter  between  profiles,  and  normalized  dissipation  is 
approximately  50%  larger  than  for  the  unstable  case  at  the  canopy  top.  In  all  cases  the 
SGS  dissipation  is  negligible  (less  than  10%  of  canopy  top  value)  below  zjh  «  0.60 . 

The  results  for  the  dense  canopy  simulations  are  given  in  Figure  5b.  The  profiles  are 
similar  to  the  sparse  canopy  with  the  dissipation  in  the  near  neutral  simulations 
approximately  50%  larger  than  the  dissipation  in  the  unstable  simulation  at  z  =  h.  One 
notable  difference  is  that  there  is  less  consistency  between  the  near  neutral  simulations 
where  SGS  dissipation  is  the  largest  for  the  high  wind  speed  simulation  (S6).  In  addition. 


the  dissipation  attenuates  more  rapidly  in  the  dense  canopy,  and  is  essentially  zero  below 
zlh  «  0^5 .  This  is  consistent  with  Shaw  and  Seginer  (1985)  who  found  the  dissipation  to 

be  negligible  within  a  dense  com  canopy  (LAI «  4.0 ),  while  dissipation  accounted  for  one- 
third  of  the  loss  in  a  sparse  artificial  wind  tunnel  canopy  (LAI «  05  ). 

Shaw  and  Seginer  (1985)  also  found  that  the  action  of  plant  canopy  drag  was  the 
primary  destmction  process  within  their  dense  canopy.  Similar  results  are  found  here 
(Figure  6).  The  general  features  show  that  the  peak  destruction  occurs  at  z//i«0.85  for 

all  simulations,  and  destmction  occurs  primarily  in  the  upper  half  of  the  canopy.  Noted 
differences  due  to  the  canopy  densities  are  larger  peak  destmction  in  the  dense  canopy 
simulations  (Figure  6b),  and  the  magnitude  decreasing  slightly  more  rapidly  within  the 
forest. 

The  primary  differences  occur  between  the  near  neutral  simulations  and  the  unstable 
simulations.  The  near  neutral  simulations  show  remarkable  agreement  throughout  the 
canopy  with  a  difference  between  profiles  of  less  than  5%.  The  unstable  simulations  show 
greater  destmction  throughout  most  of  the  canopy.  Leclerc  et  al.  (1990)  found  that 
canopy  drag  was  small  in  the  lower  portion  of  the  forest,  and  approached  zero  with 
increasing  stability  in  the  lower  half  of  the  forest. 

The  transport  terms  of  the  resolved  scale  TKE  budget  can  neither  create  nor  destroy 
energy;  rather,  the  transport  mechanisms  simply  redistribute  energy  from  one  region  to 
another.  The  general  features  of  turbulent  and  pressure  transport  have  already  been 
presented  in  Figure  2.  Here  the  details  between  simulations  are  presented,  and  the  focus 
will  be  on  the  physical  implications  of  the  profiles  near  the  canopy  top  and  below. 


The  simulated  turbulent  transport  profiles  (Figure  7)  show  a  loss  at  the  canopy  top  and 
a  gain  within  the  canopy.  The  peak  loss  for  all  simulations  occurs  at  z//i=1.05 ,  and  the 

average  loss  is  roughly  35%  of  the  combined  peak  production  at  z//i=0.95 .  Raupach  et 
al.  (1986)  found  that  turbulent  transport  loss  was  approximately  one-third  of  the 
combined  production  for  the  same  level,  Z’^h .  The  difference  here  is  that  the  simulated 
turbulent  transport  profiles  show  a  larger  gradient  in  the  upper  canopy  than  the  wind 
tunnel  plots.  This  is  probably  due  to  the  difference  in  canopy  structure,  where  in  the  wind 
tunnel  the  leaf  area  density  was  constant  with  height  in  contrast  to  the  structure  in  this 
study.  The  magnitude  of  the  peak  loss  is  very  similar  for  both  simulated  canopies  showing 
an  average  loss  of  -2.5  for  the  sparse  canopy  and  -3.0  for  the  dense  canopy;  however, 
there  appears  to  be  no  consistent  ordering  of  the  transport  terms  by  the  stability 
parameter. 

In  the  canopy  layer  the  profiles  are  roughly  ordered  in  terms  of  the  stability  parameter. 
The  sparse  canopy  simulations  (Figure  7a)  show  little  variability  between  profiles  for  the 
least  unstable  simulations  (SI  and  S2),  and  the  normalized  peak  input  of  energy  is  largest 
for  these  simulations.  Furthermore,  the  peak  normalized  turbulent  transport  source 
decreases  for  increasing  instability.  Clearly  the  peak  source  is  ordered  by  hjL .  In  the 
lower  half  of  the  canopy,  the  near  neutral  simulations  approach  a  similar  value  while 
turbulent  transport  is  larger  for  the  unstable  simulation.  The  larger  source  corresponds  to 
a  larger  gradient  of  the  TKE  flux.  Leclerc  et  al.  (1990)  found  that  the  flux  of  TKE 
increased  for  thermally  unstable  conditions,  and  normalized  budget  profiles  showed  a 


larger  turbulent  transport  source  at  the  lowest  measurement  level  in  the  forest  for  unstable 
conditions. 

The  dense  canopy  simulations  (Figure  7b)  show  similar  features.  The  peak  source  is 
ordered  by  the  stability  parameter.  The  near  neutral  simulations  collapse  to  a  single 
profile,  and  turbulent  transport  is  slightly  enhanced  in  the  lower  canopy  for  the  unstable 
simulation;  however,  there  are  several  differences  in  the  turbulent  transport  profiles  due  to 
the  change  in  canopy  density.  The  turbulent  transport  peak  is,  on  average,  larger  and 
occurs  at  a  higher  level  in  the  dense  canopy.  The  dense  canopy  profiles  attenuate  more 
rapidly  with  decreasing  height,  becoming  negligible  in  the  lowest  one-third  of  the  canopy. 
These  features  are  consistent  with  increased  normalized  shear  production  for  these 
simulations  resulting  in  greater  transport  in  the  upper  canopy,  and  the  inability  of  vertical 
motions  to  penetrate  into  the  dense  canopy. 

The  turbulent  transport  terms  plotted  in  Figures  7a  and  7b  contain  two  components. 
The  first  component  is  the  resolved  scale  turbulent  transport  ( 7] ),  and  the  second  is  the 

subgrid  scale  diffusion  of  resolved  scale  TKE.  Figure  8  gives  the  components  of  the 
turbulent  transport  for  simulation  S4,  where  7T  represents  the  sum  of  resolved  scale 
turbulent  transport  and  subgrid  scale  diffusion.  In  addition,  the  subgrid  scale  diffusion  is 
spilt  between  normal  and  shear  components  represented  as  S„  {=  -3  and 

{j=-^  respectively.  The  resolved  scale  turbulent  transport,  T^,  is  the  largest 

component  with  subgrid  scale  diffusion  generally  less  than  15%  of  the  total  (77).  Near 
the  surface  the  normal  component  of  subgrid  scale  diffusion  is  essentially  zero,  the  shear 
component  becomes  negative,  while  the  resolved  scale  turbulent  transport  is  small  and 


positive.  All  simulations  show  these  features  where  the  negative  values  of  turbulent 
transport  are  solely  due  to  the  shear  components  of  subgrid  scale  diffusion,  and  the 
negative  values  are  probably  due  to  the  subgrid  scale  parameterizations;  however,  further 
investigation  is  needed. 

An  encouraging  aspect  of  the  turbulent  transport  profiles  simulated  here  is  that  the 
results  are  well  supported  by  other  studies.  Field  experiments  for  a  pine  forest  (Lesnik, 
1974)  and  from  deciduous  forests  (Shi  et  al,  1987;  Leclerc  et  al.,  1990;  Meyers  and 
Baldocchi,  1991)  all  show  this  feature.  This  general  profile  has  also  been  found  in  wind 
tunnel  studies  (Raupach  et  al.,  1986;  Brunet  et  al,  1994)  and  in  higher  order  closure 
models  (Wilson  and  Shaw,  1977;  Meyers  and  Paw  U,  1986,  1987;  Meyers  and  Baldocchi, 
1991).  On  the  other  hand,  the  second  transport  mechanism,  pressure  transport,  is  not  as 
strongly  supported  in  the  literature. 

The  simulated  effects  of  pressure  transport  (Figure  9)  exhibit  similar  features  as 
turbulent  transport.  Pressure  transport  is  a  loss  between  0.8<z//t<13  and  a  gain  below 

z/^«0.8;  thus,  for  plant  canopy  turbulence,  pressure  transport  also  serves  to  transfer 

energy  away  from  the  primary  production  region,  the  upper  canopy.  Another  similarity 
between  transport  profiles  is  that  the  peak  loss  at  the  canopy  top  and  peak  gain  in  the 
canopy  are  larger  and  occur  at  a  higher  level  for  the  dense  canopy  simulations.  The 
pressure  transport  profiles  are  not  ordered  by  /i/L  in  the  upper  canopy,  but,  below 

z//i»0.6  (sparse  canopy)  or  z//i«0.8  (dense  canopy),  the  most  unstable  cases  show 
substantially  greater  values  of  pressure  transport. 


A  significant  difference  between  transport  turbulent  and  pressure  transport  is  that 
pressure  transport  does  not  attenuate  as  rapidly  with  decreasing  depth  into  the  canopy.  In 
addition,  pressure  transports  for  the  near  neutral  simulations  in  both  canopies  collapse  to  a 
similar  profile  below  .  In  the  lower  canopy  the  primary  balance  of  the  resolved 

scale  TBCE  budget  is  between  the  loss  due  to  canopy  drag  effects  and  the  gain  due  to 
turbulent  and  pressure  transport  with  pressure  transport  accounting  for  most  of  the  kinetic 
energy  gain.  Table  3  lists  the  average  ratio  of  resolved  scale  turbulent  transport  to 
pressure  transport  from  the  level  where  pressure  transport  first  exceeds  turbulent 
transport.  Pressure  transport  accounts  for  approximately  60%  of  the  transport  in  the 
lowest  third  of  the  sparse  canopy,  and  for  approximately  70%  of  the  transport  in  the  lower 
half  of  the  dense  canopy.  This  result  supports  the  belief  that  pressure  transport  is  not 
negligible  in  the  TKE  budget  and  that  pressure  transport  is  a  source  of  TKE  in  the  canopy 
(e.  g.  Maitani  and  Seo,  1985;  Shaw  and  Zhang,  1992;  Raupach  et  al.,  1986;  Leclerc  et  al, 
1990);  however,  this  result  contradicts  the  results  of  Brunet  et  al.  (1994)  where  pressure 
transport  was  estimated  to  be  a  loss  within  the  canopy. 

Brunet  et  al.  (1994)  calculated  the  profiles  of  shear  production  and  turbulent  transport 
from  measured  wind  tunnel  data  for  an  artificial  wheat  canopy.  The  dissipation  was  found 
using  the  measured  u  spectra,  and  a  wake  production  term  was  calculated  using  the 
formulation  presented  in  Raupach  et  al.  (1986).  Pressure  transport  was  believed  to  be  the 
dominant  term  in  the  residual.  The  residual  showed  a  gain  just  above  the  canopy  and  a 
significant  loss  within  the  canopy  indicating  a  severely  unbalanced  budget.  Brunet  et  al. 
proposed  that  a  “wake  dissipation”  term  should  be  included  in  their  original  budget  based 


on  the  assumption  that  their  measurements  could  not  resolve  the  wake  production 
process.  The  wake  dissipation  term  was  equal  and  opposite  to  the  wake  production  and  it 
was  added  to  the  viscous  dissipation,  but  their  wake  dissipation  term  only  accounted  for 
the  dissipation  of  “wake  turbulence”  by  the  mean  wind  and  not  for  the  turbulent  energy 
destroyed  by  canopy  drag.  The  new  dissipation  rate  (dissipation  =  viscous  +  wake 
dissipation)  gave  results  very  similar  to  Raupach  et  al.  (1986).  However,  the  new  residual 
still  showed  pressure  transport  as  a  TKE  sink  within  the  canopy,  and  this  new  residual  was 
almost  identical  to  the  semi-empirical  relationships  of  Zeman  and  Lumley  (1976)  and 
Deardorff  (1973).  Brunet  et  al.  (1994)  noted  that  this  consistency  cannot  be  taken  as 
proof  “since  the  closure  assumptions  are  mostly  empirical”;  however,  no  explanation  was 
given  as  to  the  observed  agreement. 

The  agreement  between  the  surface  layer  pressure  transport  parameterization  and  the 
new  residual  may  not  be  valid  for  the  plant  canopy  because  the  semi-empirical 
relationships  were  derived  for  a  smooth,  constant  flux  boundary  layer.  Both  conditions 
are  violated  in  the  plant  canopy;  however,  a  possible  test  of  the  of  the  semi-empirical 
relationships  would  be  to  calculate  the  Deardorff  (1973)  parameterization  with  the  LES 
data. 

Finally,  the  residual  term  (Figure  10)  in  the  resolved  scale  TICE  budget  represents  any 
numerical  errors,  or  an  insuflBcient  number  of  data  sets  used  to  form  the  ensemble  average 
of  the  budget  terms.  The  residual  is  largest  at  the  canopy  top  where  the  other  terms  of  the 
budget  are  large,  but  on  average  the  maximum  residual  is  less  than  2%  of  the  of  the  shear 
production  for  the  corresponding  level  in  the  upper  canopy.  This  small  residual  is 


encouraging  and  indicates  that  numerical  errors  are  negligible,  and  that  the  number  of  data 
sets  averaged  are  sufficient  to  determine  the  terms  of  the  resolved  scale  TKE  budget.  In 
addition,  the  assumptions  of  steady  state  and  horizontal  homogeneity  have  been  satisfied, 
and  it  appears  that  not  including  a  wake  production  term  in  the  SGS  kinetic  energy 
equation  has  had  little  effect  on  the  resolved  scales.  A  possible  reason  for  this  is  that 
eddies  produced  due  to  the  wake  production  process  dissipate  rapidly  (Raupach  and 
Shaw,  1982;  Wilson,  1988). 


4.  Summaiy  and  Conclusions 

LES  output  has  been  used  to  calculate  all  terms  of  the  resolved  scale  TKE  budget  for 
airflow  within  and  above  a  forest  canopy.  The  TKE  budget  contains  the  traditional  terms 
that  represent  shear  and  buoyant  production,  turbulent  transport,  and  pressure  transport; 
however,  differences  arise  because  the  budget  is  for  the  resolved  scales  of  a  LES.  These 
scales  are  mathematically  defined  using  a  numerical  filter  process  that  separates  the 
resolved  scales  of  motions  from  the  subgrid  scale  motions.  A  formulation  in  this  manner 
sets  a  lower  limit  on  the  size  of  the  resolvable  eddies.  As  a  consequence,  the  wake 
production  term  which  is  commonly  found  in  a  plant  canopy  budget  does  not  appear  in  the 
resolved  scale  budget  presented  here;  rather,  the  budget  contains  a  term  that  represents  a 
direct  extraction  of  turbulent  energy  by  the  action  of  plant  canopy  drag.  In  addition,  the 
traditional  viscous  dissipation  effects  are  represented  as  a  transfer  of  energy  from  the 
resolved  scales  to  the  subgrid  scales  by  a  parameterized  shear  stress  term. 

Several  simulations  were  run  using  different  environmental  forcing  and  plant  canopy 
specifications.  The  environmental  forcing  was  specified  through  a  mean  pressure  gradient 
that  was  adjusted  during  a  simulation  to  ensure  that  a  specified  mass  flux  was  maintained 
across  the  upwind  boundary.  The  plant  canopy  was  defined  as  either  a  sparse  canopy 
(LAI  =  2)  or  a  dense  canopy  (LAI  =  5),  and  the  effects  of  canopy  heating  on  the 
surrounding  air  were  specified  using  a  canopy  top  heat  flux.  A  canopy  heating  rate  profile 
was  calculated  fi’om  this  heat  flux  to  represent  a  high  or  a  low  solar  radiation  load  on  the 
forest.  The  thermal  stability  of  a  simulation  was  defined  by  a  stability  parameter  of  the 


form  hiL ,  where  h  is  the  canopy  top  and  L  is  the  Monin-Obukhov  length;  however, 
there  was  not  always  a  clear  order  of  the  budget  terms  by  the  stability  parameter. 

The  effect  of  canopy  density  changes  were  seen  in  all  profiles.  The  simulated  profiles 
showed  that  the  normalized  peak  value  of  a  budget  term  was  larger  in  the  dense  canopy 
simulations.  For  example,  normalized  shear  production  for  the  dense  canopy  simulations 
were  approximately  30%  larger  than  the  sparse  canopy  simulations.  This  feature  was  also 
observed  in  the  normalized  profiles  of  buoyant  production,  canopy  drag  effects, 
dissipation,  and  the  transport  terms.  In  addition,  the  dense  canopy  simulations  showed 
that  the  magnitude  of  the  budget  terms  attenuated  much  more  rapidly  within  the  canopy. 
This  indicated  that  turbulent  motions  are  suppressed  in  the  lower  portions  of  a  dense 
forest  when  compared  to  a  more  open  forest. 

The  most  significant  aspect  of  this  study  was  that  the  three  dimensional  pressure 
perturbation  field  was  numerically  calculated;  thus,  allowing  the  calculation  of  the  pressure 
transport  term  at  all  levels  in  the  domain,  and  the  effects  of  pressure  transport  could  now 
be  compared  to  the  other  budget  terms  for  the  first  time.  The  general  features  of  the 
budget  showed,  that  above  the  canopy,  there  was  an  approximate  balance  between 
production  and  subgrid  scale  transfer  (dissipation),  and  there  were  indications  that 
turbulent  and  pressure  transport  roughly  opposed  each  other  with  pressure  transport 
balancing  approximately  30%  of  turbulent  transport;  however,  there  was  considerable 
variation  between  simulations.  This  approximate  balance  did  not  exist  near  the  canopy  top 
and  below.  At  the  canopy  top,  production  was  balanced  by  plant  canopy  drag  effects,  the 
energy  transfer  to  the  subgrid  scales  (dissipation),  and  by  both  transport  mechanisms. 


where  turbulent  and  pressure  transport  counteracted  over  half  of  the  production  in  the 
simulations.  In  the  lowest  two-thirds  of  the  canopy,  production  and  subgrid  scale  transfer 
became  negligible,  and  the  primary  balance  was  between  the  loss  due  to  canopy  drag 
effects  and  the  gain  due  to  transport.  Further,  pressure  transport  was  the  dominant  source 
in  the  lowest  levels  of  the  forest.  In  fact,  pressure  transport  accounted  for  approximately 
60%  of  the  TKE  source  in  the  sparse  canopy  and  roughly  70%  of  the  source  in  the  dense 


canopy. 


Table  1.  Plant  canopy  variables  and  forcing  for  simulations. 


Simulation  LAI  Q*  U  hjL  u»  Data 

(ms‘‘  K)  (ms'*)  (ms'*)  Sets 


SI 

2.0 

0.005 

2.0 

-0.024 

0.388 

15 

S2 

2.0 

0.125 

4.0 

-0.074 

0.767 

13 

S3 

2.0 

0.005 

1.0 

-0.185 

0.196 

15 

S4 

2.0 

0.125 

1.0 

-1.398 

0.281 

12 

S5 

5.0 

0.005 

2.0 

-0.017 

0.424 

12 

S6 

5.0 

0.125 

4.0 

-0.064 

0.801 

12 

S7 

5.0 

0.005 

1.0 

-0.153 

0.209 

06 

S8 

5.0 

0.125 

1.0 

-1.252 

0.302 

08 

Q*  =  canopy  top  heat  flux;  U  =  mean  wind  speed;  h  =  canopy  height; 

L  =  canopy  top  Monin-Obukov  length;  u*  =  canopy  top  horizontal  mean  friction  velocity 
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Table  2.  Ratio  of  turbulent  transport  ( 77^  )  and  pressure  transport  ( T  ^ ) 

peak  loss  at  the  canopy  top  to  the  combined  shear  and  buoyant  peak 
production  (  P  +  ). 


Simulation  L 


SI 

0.39 

0.13 

S2 

0.32 

0.15 

S3 

0.43 

0.11 

S4 

0.42 

0.12 

S5 

0.38 

0.13 

S6 

0.35 

0.13 

S7 

0.37 

0.14 

S8 

0.40 

0.17 
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Table  3.  Average  ratio  of  resolved  scale  turbulent  transport  (  7] )  to  pressure 
transport  ( T  ).  The  average  is  calculated  over  the  region  from  the 
surface  to  the  level  (z//i)  where  pressure  transport  first  becomes 
larger  than  resolved  scale  turbulent  transport. 


Simulation  zjh 


SI 

0.35 

0.42 

S2 

0.35 

0.45 

S3 

0.35 

0.29 

S4 

0.45 

0.43 

S5 

0.55 

0.27 

S6 

0.55 

0.30 

S7 

0.65 

0.30 

S8 

0.65 

0.25 

LAD 


Figure  1.  Vertical  distribution  of  leaf  area  density  Integration  gives  leaf  area 
index  (LAI)  of  2  and  5. 


Normalized  TKE  Budget 
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Figure  2.  Normalized  resolved  scale  TKE  budget  for  simulation  S5,  /i/L  =  -0.017 .  Plot 
labels  refer  to  budget  terms  in  equation  (13),  and  R  represents  the  residual. 


Shear  Production 


Figure  3a.  Normalized  shear  production  (  /^ )  for  a  sparse  canopy  LAI  =  2.  Plot  labels 
are  ordered  in  terms  of  increasingly  negative  hjL ;  S 1  ( h/L = -0.024 ), 

S2  ( -0.074  ),  S3  ( -0.185  ),  S4  ( -1398 ). 


Figure  3b.  Normalized  shear  production  (  Pj )  for  a  dense  canopy  LAI  =  5.  Plot  labels 
are  ordered  in  terms  of  increasingly  negative  h/L ;  S5(/i/L= -0.017), 

S6( -0.064),  S7( -0.153),  S8( -1.252). 


Buoyant  Production 


Figure  4a.  Normalized  buoyant  production  ( )  for  a  sparse  canopy  LAI  =  2.  Plot  labels 
are  ordered  in  terms  of  increasingly  negative  h/L ;  SI  (/2/L= -0.024), 

S2  ( -0.074  ),  S3  ( -0.185 ),  S4  ( -1398  ). 


Buoyant  Production 


Figure  4b.  Normalized  buoyant  production  ( )  for  a  dense  canopy  LAI  =  5.  Plot  labels 
are  ordered  in  terms  of  increasingly  negative  hjL ;  S5(^/L=-0.017), 

S6( -0.064),  S7( -0.153),  S8( -1.252). 


Figure  5a.  Normalized  SGS  dissipation  (D^gs )  for  a  sparse  canopy  LAI  =  2.  Plot  labels 

are  ordered  in  terms  of  increasingly  negative  hjL ;  S 1  ( h/L  -  -0.024  ), 

S2  (-0.074),  S3  (-0.185),  S4  (-1398). 


Figure  5b.  Normalized  SGS  dissipation  (D^gs )  for  a  dense  canopy  LAI  =  5.  Plot  labels 
are  ordered  in  terms  of  increasingly  negative  hjL ;  S5(/j/L= -0.017), 

S6( -0.064),  S7( -0.153),  S8(-U52). 


Figure  6a.  Normalized  canopy  drag  effects  {D^d  )  for  a  sparse  canopy  LAI  =  2.  Plot 

labels  are  ordered  in  terms  of  increasingly  negative  hjL ;  SI  (/i/L= -0.024), 
S2  ( -0.074 ),  S3  ( -0.185 ),  S4  ( -1398 ). 


Figure  6b.  Normalized  canopy  drag  effects  ( )  for  a  dense  canopy  LAI  =  5.  Plot 

labels  are  ordered  in  terms  of  increasingly  negative  hjL ;  S5(/i/L=-0.017 ), 
S6( -0.064),  S7( -0.153),  S8(-1252). 


Figure  8.  Normalized  components  of  the  turbulent  transport  term  ( 7T)  for  simulation 
S4.  T^  is  the  resolved  scale  component.  S„  and  5,  represent  the  normal  and  shear 
components  of  subgrid  scale  dffiision  of  resolved  scale  TKE,  respectively. 


Figure  9a.  Normalized  pressure  transport  ( 7^  )  for  a  sparse  canopy  LAI  =  2.  Plot  labels 
are  ordered  in  terms  of  increasingly  negative  hjL ;  S 1  ( h/L = -0.024  ), 

S2  (-0.074),  S3  (-0.185),  S4  (-1398). 


Figure  10a.  Normalized  residual  ( /? )  for  a  sparse  canopy  LAI  =  2.  Plot  labels  are 
ordered  in  terms  of  increasingly  negative  hjL ;  SI  (/i/L= -0.024), 

S2 (-0.074),  S3  (-0.185),  S4(-U98). 


Figure  10b.  Normalized  residual  ( /? )  for  a  dense  canopy  LAI  =  5 .  Plot  labels  are 

ordered  in  terms  of  increasingly  negative  hjL ;  S5(hlL= -0.017  ),  S6(  -0.064  ), 
S7( -0.153),  S8( -1.252). 


Appendix  A 


Derivation  of  the  resolved  scale  momentum  equation 

The  resolved  scale  momentum  equation  was  derived  from  the  Navier-Stokes  equation 
for  an  incompressible,  Boussinesq  fluid.  In  tensor  notation  the  momentum  equation  is 


dt 


dx. 


p  dx^  e. 


(Al) 


In  LES  the  resolved  scale  part  of  the  flow  is  numerically  defined  by  applying  a  filter 
function  to  the  governing  equations.  Filtering  the  Navier-Stokes  equation  gives 


du.  ^  ^u.u.  \  dp  ^  g 

dt  dx.  pdx^  9^ 


+  —BS.,  +  V- 


d^u. 


dx.dXj 


(A2) 


where  the  overbar  represents  the  filter  process.  The  velocity  components  in  the  advection 
term  must  be  partitioned  into  a  resolved  scale  part  and  a  subgrid  scale  part,  u.=u.+  u'. , 


to  avoid  writing  an  equation  for  u.u.  (Leonard,  1974).  In  addition,  the  presence  of  the 

plant  canopy  complicates  the  filtering  process.  The  horizontal  pressure  differentiation  and 
the  Laplacian  should  not  commute  with  the  filtering  due  to  discontinuities  across  the  plant 
elements  according  to  arguments  concerning  horizontal  averaging  and  differentiation 
presented  in  Wilson  and  Shaw  (1977)  and  Raupach  and  Shaw  (1982).  Therefore,  the 
pressure  was  also  partitioned  into  a  resolved  scale  part  and  a  subgrid  scale  part,  and 
performing  the  substitution  gives 
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dUM.  ^(u;Uj+UM'.+ulu>.) 


dt  dx. 


dx. 

j 


p  dx.  p  dx. 


(A3) 


g  d^U. 

+■^95..,  +  y - !- 

6>,  ■’  dx.dx. 


d^u' 


+  V- 


dx.dx. 

/  t 


where  the  second  term  on  the  right  hand  side  represents  the  energy  transfer  between  the 
resolved  and  the  subgrid  scales.  The  fourth  and  seventh  terms  represent  form  and  viscous 
drag  across  the  plant  elements,  respectively  (Wilson  and  Shaw,  1977;  Raupach  and  Shaw, 
1982).  Both  effects  were  parameterized  as  one  term,  and  for  ease  of  presentation  the 
viscous  drag  term  will  not  be  carried  through  in  the  derivation.  In  addition,  the  effects  of 
molecular  diffiision  are  negligible  on  the  resolved  scale  and  are  not  retained. 

The  filtering  process  is  different  from  the  traditional  statistical  averaging  for  turbulent 
flow.  In  general. 


/'  0  and 


(A4) 


However,  the  inequalities  in  (A4)  could  be  equated  depending  on  the  chosen  filter 
(Reynolds,  1990).  The  filter  used  in  the  LES  here  was  chosen  such  that  the  equalities 
hold;  thus,  (A3)  becomes 


dt  dx.  dx.  p  dx.^  9„  P  dx.^ 


(A5) 


Next,  the  advection  term  in  current  LES  were  written  in  rotational  form  to  give 
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Equation  (A7)  represents  the  subgrid  scale  shear  stresses,  where  the  normal  components 
were  subtracted,  and  an  equal  and  opposite  term  was  added  to  the  resolved  scale  pressure 
(A6).  In  addition,  the  horizontal  mean  pressure  (angled  brackets)  gradient  was  separated 
from  the  resolved  scale  pressure  and  treated  as  an  external  forcing.  As  a  result  of  the 
separation,  the  pressure  in  the  fourth  term  on  the  right  hand  side  represents  the  deviation 
from  the  horizontal  mean  indicated  by  the  double  primes.  The  resolved  scale  momentum 
equation  that  was  used  to  form  the  turbulent  kinetic  energy  budget  is 


II 

1 

^  du.  du^ 
1  1 

d 

fu.u.'\ 

J  J 

dt  ^ 

.dx.  dx.j 

dx^ 

^  2  > 

2  g  nc 


(A9) 


1  <7  fp"^ 


pdx. 


1^(P)  \dp> 


\PJ 


p  dx.  dx.  pdx. 

where  e'  =  u'm'Ji  is  the  subgrid  scale  kinetic  energy. 

Finally,  two  parameterizations  were  used  in  (A9).  The  subgrid  scale  shear  stress  was 
represented  in  terms  of  the  resolved  scale  strain  rate  of  the  form 
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where  the  difiEusion  coefficient  was  represented  as  a  function  of  the  subgrid  scale  kinetic 
energy  and  a  dissipation  length  (Moeng,  1984). 

The  canopy  drag  effects  were  parameterized  as 


pdx. 


(All) 


where  Cd  is  the  isotropic  drag  coefficient;  a(z)  is  the  leaf  area  density;  and  V  is  the  vector 


magnitude  wind. 
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Appendix  B 

Derivation  of  the  horizontal  mean  momentum  equation 

The  horizontal  mean  momentum  equation  was  derived  from  the  resolved  scale 
momentum  equation  (A9)  in  the  following  manner:  first  the  flow  variables  were  portioned 
into  a  horizontal  mean  and  a  deviation  therefrom;  such  as,  «,  =(«,)  +  u" ;  and  a 

horizontal  average  was  then  performed  to  give  the  horizontal  mean  momentum  equation. 
A  term  by  term  derivation  follows: 


The  left  hand  side  of  (A9)  is 


where  the  time  differentiation  and  horizontal  average  commute  provided  the  plant  parts 
are  not  waving  (Raupach  and  Shaw,  1982),  =  0  by  definition.  On  the  right  hand  side 


of  (A9)  the  expansion  of  the  rotational  part  of  the  advection  term  is 
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The  averaging  process  and  horizontal  differentiation  commute  in  this  case  because  the 
first-order  spatial  differentiation  and  horizontal  averaging  commute  (Raupach  and  Shaw, 
1982).  As  a  result,  the  second  and  third  term  are  zero  on  the  horizontal  average. 
Partitioning  the  kinetic  energy  part  of  the  advection  term  gives 
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where  the  second  term  drops  out  on  the  horizontal  average.  The  buoyancy  term  is 


The  plant  canopy  effects,  form  drag,  are  parameterized  in  the  model  in  the  following 
manner 

--^  =  -C^aiz)m  (B5) 

I 

and  this  is  symbolically  partitioned  as 
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The  SGS  kinetic  energy  gradient  is  partitioned  in  a  similar  manner,  and  it  equals 
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Finally  partitioning  the  SGS  shear  stresses  and  performing  the  horizontal  average  gives 


The  eddy  difiusivity  is  a  fimction  of  the  SGS  kinetic  energy;  thus,  the  horizontal  averaging 
and  differentiation  should  commute.  Therefore,  the  second  and  third  terms  will  equal 
zero,  and  the  result  is 


Finally,  combining  terms  gives  the  resolved  scale  horizontal  mean  momentum 


vni 


where  the  terms  represent 

Term  I:  local  time  rate  of  change  of  horizontal  mean  momentum 

Term  II:  advection  of  horizontal  mean  momentum  by  the  mean  wind 

Term  III:  Reynolds  stress  terms 

T  erm  IV :  buoyancy 

Term  V:  mean  pressure  gradient 

Term  VI:  form  drag  on  the  mean  wind  due  to  plant  canopy-airflow  interaction 

Term  VII:  SGS  normal  stresses 

Term  VIII:  parameterized  SGS  shear  stresses 


Appendix  C 


Derivation  of  the  horizontal  deviation  momentum  equation 

The  horizontal  deviation  momentum  equation  was  derived  in  the  conventional  manner 
where  the  horizontal  mean  momentum  equation  (BIO)  was  subtracted  from  the  resolved 
scale  momentum  equation  (A9).  The  flow  variables  in  the  resolved  scale  momentum 
equation  were  first  partitioned  into  a  horizontal  mean  and  a  deviation  therefrom  of  the 
form: 


and  subtract  the  horizontal  mean  momentum  equation  giving 
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The  advection  terms  can  be  simplified  using  the  product  rule,  and  equation  (C2)  can  be 
written  as 


where  the  SGS  stresses  are  defined  as 
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The  terms  of  (C3)  represent 

Term  I:  local  time  rate  of  change  of  deviation  velocity  component 

Term  II:  advection  of  the  deviation  velocity  by  the  mean  wind 
Term  III:  mean  shear  term 

Term  IV :  advection  of  the  deviation  velocity  by  deviation  velocity 

T  erm  V :  Reynolds  stress  term 

Term  VI:  buoyancy 

Term  VII:  deviation  pressure  gradient 

Term  VIII:  plant  canopy  drag  effects 

Term  IX:  SGS  normal  stresses 

Term  X:  parameterized  SGS  shear  stress 
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Appendix  D 

Derivation  of  the  resolved  scale  turbulent  kinetic  energy  equation 

The  resolved  scale  TKE  equation  for  canopy  flow  was  derived  by  multiplying  the 
horizontal  deviation  momentum  equation  (C3)  by  a  horizontal  deviation  velocity 
component,  u" ,  and  then  performing  a  horizontal  average.  The  result  is: 
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For  steady,  horizontally  homogeneous  conditions  with  negligible  mean  subsidence,  the 
TKE  budget  becomes 
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where  E"  =  u!'u!'ll  is  the  resolved  scale  turbulent  kinetic  energy,  and  the  terms  represent 
the  following:  P  is  the  mean  shear  production;  is  the  buoyant  production;  7]  is  the 
turbulent  transport;  7  is  the  pressure  transport;  is  the  parameterized  influence  of 
plant  canopy  drag;  and  represents  the  combined  effects  of  subgrid  scale  diffusion  of 
resolved  scale  TKE  (first  two  terms)  and  the  transfer  of  resolved  scale  TKE  to  the  subgrid 


scales  (last  term). 
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